#include <iostream>
#include <fstream>
#include <cstring>
#include <ctime>
#include <cstdlib>
#include "ti.h"
#include "tplate.h"
#include "txtpara.h"
#include "regression.h"
using namespace std;

int main(int argc, char** argv)
{
    if (argc != 5)
    {
        cout << "Usage:\n" << argv[0] <<
            " param_file training_image_file template_file paramfile"
            << endl;
        return 0;
    }
    float* ti;
    int length, width, NR_variable;
    read_image_cont(&ti, argv[2], length, width, NR_variable);
    float rate = 1;
    rate = getparaFLOAT("rate", argv[1]);
    tplate tp;
    // Read the training image to the trianing grid
    load_template_3D(&tp, argv[3]);
    //Randomly select 10% of the points
    srand(time(NULL));
    int nr_sample = length * width * rate;
    float* data = new float[4 * tp.tpsize * nr_sample * NR_variable];
    float* ydata = new float[tp.tpsize * nr_sample * NR_variable];
    float* theta = new float[8 * tp.tpsize * NR_variable];    //w0, w1, w2, w3, and sigma, and pvalues for w1, w2, w3
    for (int k = 0; k < NR_variable; k++)
        for (int i = 0; i < tp.tpsize; i++)
            for (int j = 0; j < nr_sample * 4; j++)
                data[k * 4 * tp.tpsize * nr_sample + i * nr_sample * 4 + j] = 0;
    //The first column remains 0
    for (int i = 0; i < nr_sample; i++)
    {
        int xpos = rand() % length;
        int ypos = rand() % width;
        bool out_bound = false;
        for (int tpos = 0; tpos < tp.tpsize; tpos++)
        {
            int x_ = xpos + tp.pos[2 * tpos];
            int y_ = ypos + tp.pos[2 * tpos + 1];
            if ((x_ < 0) || (x_ >= length))
            {
                out_bound = true;
                break;
            }
            if ((y_ < 0) || (y_ >= width))
            {
                out_bound = true;
                break;
            }
            data[tpos * nr_sample * 4 + i * 4 + 0] = 1;
            data[tpos * nr_sample * 4 + i * 4 + 1] = ti[x_ * width + y_];
            data[tpos * nr_sample * 4 + i * 4 + 2] = ti[x_ * width + y_ + width * length];
            data[tpos * nr_sample * 4 + i * 4 + 3] = ti[x_ * width + y_ + width * length * 2];
            ydata[nr_sample * tpos + i] = ti[xpos * width + ypos];

            data[4 * tp.tpsize * nr_sample + tpos * nr_sample * 4 + i * 4 + 0] = 1;
            data[4 * tp.tpsize * nr_sample + tpos * nr_sample * 4 + i * 4 + 1] = ti[x_ * width + y_];
            data[4 * tp.tpsize * nr_sample + tpos * nr_sample * 4 + i * 4 + 2] = ti[x_ * width + y_ + width * length];
            data[4 * tp.tpsize * nr_sample + tpos * nr_sample * 4 + i * 4 + 3] = ti[x_ * width + y_ + width * length * 2];
            ydata[tp.tpsize * nr_sample + nr_sample * tpos + i] = ti[xpos * width + ypos + width * length];

            data[4 * tp.tpsize * nr_sample * 2 + tpos * nr_sample * 4 + i * 4 + 0] = 1;
            data[4 * tp.tpsize * nr_sample * 2 + tpos * nr_sample * 4 + i * 4 + 1] = ti[x_ * width + y_];
            data[4 * tp.tpsize * nr_sample * 2 + tpos * nr_sample * 4 + i * 4 + 2] = ti[x_ * width + y_ + width * length];
            data[4 * tp.tpsize * nr_sample * 2 + tpos * nr_sample * 4 + i * 4 + 3] = ti[x_ * width + y_ + width * length * 2];
            ydata[tp.tpsize * nr_sample * 2 + nr_sample * tpos + i] = ti[xpos * width + ypos + width * length * 2];
        }
        if (out_bound == true)
            i -= 1;
    }
    ofstream param;
    param.open(argv[4]);
    param << 8 << " ";            //the size of each theta
    param << tp.tpsize << endl;    //the total number os thetas
    for (int tpos = 0; tpos < tp.tpsize; tpos++)
        for (int k = 0; k < NR_variable; k++)
        {
            //Calculate regression coefficient
            regression(data + 4 * nr_sample * tpos + 4 * tp.tpsize * nr_sample * k, 4, nr_sample,
                ydata + nr_sample * tpos + tp.tpsize * nr_sample * k, 1, theta + 8 * tpos + 8 * tp.tpsize * k);
            for (int i = 0; i < 8; i++)
            {
                param << theta[8 * tp.tpsize * k + 8 * tpos + i] << " ";    // the content of each theta
            }
            param << endl;
        }
    param.close();
    return 0;
}
